﻿// mod09ga_prepare.cpp : 此文件包含 "main" 函数。程序执行将在此处开始并结束。
//使用C盘 OSGeo4W的GDAL库快速开发GDAL程序。
// 用在wyj老师的气候中心风三D自研NDVI合成工作里。
// 这个程序用来判识云下阴影，对于云的判断主要根据蓝通道阈值，这个蓝通道可以TOA也可以是SR，反正阈值需要作为输入参数。
// 输入文件必须是WGS84 坐标。
// 1.这里面第一步用蓝通道阈值判断云，这里没有考虑积雪什么的，反正过亮的像素一律作为云。
// 2.然后把1的云像素为中心扩展3x3 box都假设为云。
// 3.针对2中的云，假设云高2km，尽管这个数字肯定不准确，但是前面也说了3x3都假设是云，也就是说基本能把云高的误差减到最低了。
// 4.使用T1、P1计算云的真实位置RC。
// 5.使用T0、P0和RC计算阴影的位置SH。
// 6.输出结果 0未处理、1云、2云阴影。
// 2024/5/12
// 
// This is for linux 
// 2024/5/20

#include <memory>
#include <iostream>
#include "gdal_priv.h"
#include <vector>
#include <string>
#include "getopt_pp.h"
#define WGDALRASTER_H_IMPLEMENTATION
#include "../sharedcodes/wGdalRaster.h"
using namespace GetOpt;
using namespace std;

const float ang2rad = 3.14 / 180.0;
const int NONE = 0;
const int CLOUD = 1;
const int SHADOW = 2;

struct Posi2D {
    Posi2D() :x(0), y(0) {}
    Posi2D(int xx, int yy) :x(xx), y(yy) {}
    int x, y;
};

void getRealCloudRelativeXY(float viewAzimuthAngle, float dist, float& dx,float& dy)
{
    float p2 = viewAzimuthAngle;
    if (viewAzimuthAngle <= 90.0) {
        dx = sinf(p2 * ang2rad) * dist;
        dy = -cosf(p2 * ang2rad) * dist;
    }
    else if (viewAzimuthAngle <= 180.0) {
        p2 = 180.f - viewAzimuthAngle;
        dx = sinf(p2 * ang2rad) * dist;
        dy = cosf(p2 * ang2rad) * dist;
    }
    else if (viewAzimuthAngle <= 270.0) {
        p2 = viewAzimuthAngle - 180.f;
        dx = -sinf(p2 * ang2rad) * dist;
        dy = cosf(p2 * ang2rad) * dist;
    }
    else {
        p2 = 360.f - viewAzimuthAngle;
        dx = -sinf(p2 * ang2rad) * dist;
        dy = -cosf(p2 * ang2rad) * dist;
    }
}

void computeShadowRelativePosition(const float cloudHeight,float solarZenithAngle, float solarAzimuthAngle ,float& dist, float& dx,float& dy) {
    dist = tanf(solarZenithAngle * 3.14 / 180) * cloudHeight;
    if (solarAzimuthAngle <= 90.0) {
        float p2 = solarAzimuthAngle;
        dx = -sinf(p2 * ang2rad) * dist;
        dy = cosf(p2 * ang2rad) * dist;
    }
    else if (solarAzimuthAngle <= 180.0) {
        float p2 = 180.0 - solarAzimuthAngle;
        dx = -sinf(p2 * ang2rad) * dist;
        dy =  -cosf(p2 * ang2rad) * dist;
    }
    else if (solarAzimuthAngle <= 270.0) {
        float p2 =  solarAzimuthAngle - 180.0 ;
        dx = sinf(p2 * ang2rad) * dist;
        dy = -cosf(p2 * ang2rad) * dist;
    }
    else {
        float p2 = 360.0 - solarAzimuthAngle;
        dx = sinf(p2 * ang2rad) * dist;
        dy = cosf(p2 * ang2rad) * dist;
    }
}


int main(int argc,char* argv[])
{
    cout << "A program to detect cloud shadows." << endl;
    cout<<"version v1.0.0.2 in linux"<<endl;
    cout << "usage:cloud_shadows  " << endl;
    cout<<"  -i inputfilename must be wgs84."<<endl;
    cout<<"  -o outputfilename value:0-unknown, 1-cloud, 2-shadow."<<endl;
    cout<<"  --blue bandindex of blue, 0 based."<<endl;
    cout<<"  --t0 bandindex of blue, 0 based."<<endl;
    cout<<"  --t1 bandindex of blue, 0 based."<<endl;
    cout<<"  --p0 bandindex of blue, 0 based."<<endl;
    cout<<"  --p1 bandindex of blue, 0 based."<<endl;
    cout<<"  --gs geo data scale."<<endl;
    cout<<"  --bt blue threshold for cloud (blue gt bt is cloud)"<<endl;
    GDALAllRegister();

    GetOpt_pp ops(argc, argv);
    


    string infilename  ;
    string outname ;
    int blueindex = -9999;
    int t0index = -9999;
    int t1index = -9999;
    int p0index = -9999;
    int p1index = -9999;
    float geoScale = -9999;
    float blueThreshold = -9999;
    

    ops >>Option('i',infilename)
        >>Option('o', outname) 
        >>Option("blue",blueindex)   
        >>Option("t0",t0index) 
        >>Option("t1",t1index) 
        >>Option("p0",p0index) 
        >>Option("p1",p1index) 
        >>Option("gs",geoScale) 
        >>Option("bt",blueThreshold)   
        ;



    cout << "input:" << endl;
    cout << "infilename "<< infilename <<endl;
    cout << "blueindex " << blueindex << endl;
    cout << "t0index " << t0index << endl;
    cout << "t1index " << t1index << endl;
    cout << "p0index " << p0index << endl;
    cout << "p1index " << p1index << endl;
    cout << "geoScale " << geoScale << endl;
    cout << "blueThreshold " << blueThreshold << endl;
    cout << "outname " << outname << endl;

    if( infilename=="") {
        cout<<"bad infilename "<<endl ;
        return 11 ;
    }
        if( outname=="") {
        cout<<"bad outname "<<endl ;
        return 11 ;
    }
    if( blueindex<0 || t0index<0 || t1index<0 || p0index<0 || p1index<0 || geoScale<0 || blueThreshold<0 ) {
        cout<<"some values invalid."<<endl;
        return 11 ;
    }

    shared_ptr<wGdalRaster> inRaster( wGdalRasterFactory::OpenFile(infilename) );
    if (inRaster.get() == 0) {
        cout << "failed to open " << infilename;
        return 12;
    }

    int xsize = inRaster->getXSize();
    int ysize = inRaster->getYSize();
    int asize = xsize * ysize;
    
    wGdalRasterByte outCloudShadow;
    outCloudShadow.create(xsize, ysize, 1);
    outCloudShadow.fill(0, NONE );
    outCloudShadow.copyProj(inRaster->getProj());
    outCloudShadow.copyTrans(inRaster->getTrans());

    float cellSizeInKm = 1.0;
    float cloudHeight = 2.0;//in km.
    int per1 = -1;
    //第一轮筛选云
    for (int it = 0; it < asize; ++it) {
        float tblue = inRaster->getValuef(it, blueindex);
        if (tblue > 800 && tblue < 20000) {
            outCloudShadow.setValuei(it, 0, CLOUD);
            int iy = it / xsize;
            int ix = it % xsize;
            for (int dx = -2; dx <= 2; ++dx) {
                for (int dy = -2; dy <= 2; ++dy) {
                    int it2 = (iy + dy) * xsize + ix + dx;
                    float tblue2 = inRaster->getValuef(it2, blueindex);
                    if (tblue2 > 400 && tblue2 < 20000) {
                        if (outCloudShadow.getValuei(it2, 0) == NONE) {
                            outCloudShadow.setValuei(it2, 0, CLOUD);
                        }
                    }
                }
            }
        }
    }

    //第二轮检查阴影
    for (int it = 0; it < asize; ++it) {
        float tblue = inRaster->getValuef(it, blueindex);
        float tt0 = inRaster->getValuef(it, t0index) * 0.01f ;
        float tt1 = inRaster->getValuef(it, t1index) *0.01f;
        float tp0 = inRaster->getValuef(it, p0index)*0.01f;
        float tp1 = inRaster->getValuef(it, p1index)*0.01f;
        int firstOutValue = outCloudShadow.getValuei(it, 0);
        int iy = it / xsize;
        int ix = it % xsize;

        int outval = 0;
        if (tt0 >= 0 && tt0 < 89 && tt1 >= 0 && tt1 < 89 && tp0 >= 0 && tp0 < 360.1 && tp1 >= 0 && tp1 < 360.1) {
            if (firstOutValue == CLOUD) {
                {
                    // 这个800的阈值表示较亮的云，应该在这个像素周围再扩2个像素作为潜在云边缘，然后参与阴影检验。
                    //is cloud
                    //try to guess it's shadow
                    float distCloudPixel2CloudReal = cloudHeight * tanf(tt1 * ang2rad);
                    float realCloudDx = 0;
                    float realCloudDy = 0;
                    getRealCloudRelativeXY(tp1, distCloudPixel2CloudReal, realCloudDx, realCloudDy);
                    int realCloudX = ix + realCloudDx - 0.5 ;
                    int realCloudY = iy + realCloudDy - 0.5 ;
                    //make a box of 2x2 for possible real cloud positions.
                    vector<Posi2D> possRealCloudPosiArr;
                    possRealCloudPosiArr.push_back(Posi2D(realCloudX, realCloudY));
                    possRealCloudPosiArr.push_back(Posi2D(realCloudX + 1, realCloudY));
                    possRealCloudPosiArr.push_back(Posi2D(realCloudX + 1, realCloudY + 1));
                    possRealCloudPosiArr.push_back(Posi2D(realCloudX, realCloudY + 1));

                    float tdist2Shadow = 0;
                    float shadowDx = 0;
                    float shadowDy = 0;
                    computeShadowRelativePosition(cloudHeight, tt0, tp0, tdist2Shadow, shadowDx, shadowDy);
                    vector<Posi2D> possShadowPosiArr;
                    for (int ic = 0; ic < possRealCloudPosiArr.size(); ++ic) {
                        int shadow_x = possRealCloudPosiArr[ic].x + shadowDx - 0.5;
                        int shadow_y = possRealCloudPosiArr[ic].y + shadowDy - 0.5;

                        possShadowPosiArr.push_back(Posi2D(shadow_x,shadow_y));
                        possShadowPosiArr.push_back(Posi2D(shadow_x + 1, shadow_y));
                        possShadowPosiArr.push_back(Posi2D(shadow_x + 1, shadow_y + 1));
                        possShadowPosiArr.push_back(Posi2D(shadow_x, shadow_y + 1));
                    }
                    for (int ish = 0; ish < possShadowPosiArr.size(); ++ish) {
                        int shx = possShadowPosiArr[ish].x;
                        int shy = possShadowPosiArr[ish].y;
                        if (shx >= 0 && shx < xsize && shy >= 0 && shy < ysize) {
                            if (outCloudShadow.getValuei(shx, shy, 0) == NONE ) {
                                outCloudShadow.setValuei(shx, shy, 0, SHADOW);//set shadow value
                            }
                        }
                    }

                }
            }
        }
        int per = it * 100.0 / asize;
        if (per != per1) {
            per1 = per;
            cout << per1 << endl; 
        }
    }
    cout << "100" << endl;
    bool saveok = outCloudShadow.save(outname);
    if (saveok == false) {
        cout << "failed to save " << outname << endl;
        return 13;
    }
    cout << "output in " << outname << endl;
    cout << "done" << endl;
    return 0;
}


